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A  group  of  algorithms  is  presented  generalizing  the  Fast  Fourier  Transform  to  the  case  of  non¬ 
integer  frequencies  and  non-equispaced  nodes  on  the  interval  J— The  schemes  of  this  paper 
are  based  on  a  combination  of  certain  analytical  considerations  with  the  classical  Fast  Fourier 
Transform,  and  generalize  both  the  forward  and  backward  FFTs.  Each  of  the  algorithms  requires 
0{N  -logN  +  N  •log(l/£))  arithmetic  operations,  where  £  is  the  precision  of  computations  and  N 
is  the  number  of  nodes.  The  efficiency  of  the  approach  is  illustrated  by  several  numerical  examples. 
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1  Introduction 


Fourier  techniques  have  been  a  popular  analytical  tool  in  the  study  of  physics  and  engineering 
for  more  than  two  centuries.  A  reason  for  the  usefulness  of  such  techniques  is  that  the  trigono¬ 
metric  functions  e*"®  axe  eigenfunctions  of  the  differentiation  operator  and  can  be  effectively 
used  to  model  solutions  of  differential  equations  which  arise  in  the  fields  mentioned  above. 

More  recently,  the  arrival  of  digital  computers  and  the  development  of  the  Fast  Fourier 
Transform  (FFT)  algorithm  in  the  1960s  have  established  Fourier  analysis  as  a  powerful  and 
practical  numerical  tool.  The  FFT,  which  computes  discrete  Fourier  transforms  (DFTs),  is 
now  central  to  many  areas,  most  notably  spectral  analysis  and  signal  processing.  In  some 
applications  however,  the  input  data  is  not  uniformly  spaced,  a  condition  which  is  required 
for  the  FFT.  In  this  paper  we  present  a  set  of  algorithms  for  computing  more  efficiently  some 
generalizations  of  the  DFT,  namely  the  forward  and  inverse  transformations  described  by  the 
equations 

N 

ifc=0 

for  j  =  0, . . . ,  N,  where  fj  €  C,  a*  €  C,  wjt  €  [-N/2,  N/2]  and  Xj  6  [-tt,  tt].  Each  algorithm 
requires  a  number  of  arithmetic  operations  proportional  to 

JV-logAT  +  AT-logQ)  (2) 

where  e  is  the  desired  accuracy,  compared  with  0{N^)  operations  required  for  the  direct  ap¬ 
plication  and  O(iV^)  for  the  direct  inversion. 

Remark  1.1  The  DFT  in  “unaliased”  form,  as  described  by  the  equations 

N/2-1 

u=  Y.  (3) 

k=-N/2 

is  clearly  a  special  case  of  (1);  the  FFT  algorithm  reduces  the  number  of  operations  for  the 
DFT  from  to  0{N  -log  N)  by  a  sequence  of  algebraic  manipulations.  In  the  more  general 

case  of  (1),  the  structure  of  the  linear  transformation  is  also  exploitable,  and  the  algorithms  of 
this  paper  combine  certain  analytical  results  with  the  existing  FFT. 

The  plan  of  the  paper  is  as  follows.  We  start  in  Section  2  with  some  results  from  analysis  and 
approximation  theory  which  are  used  in  the  design  of  the  algorithms.  An  exact  statement  of  the 
problem  in  Section  3  is  then  followed  by  informal  descriptions  of  the  algorithms  in  Section  4. 
In  Section  5  we  introduce  some  notation  which  is  used  in  a  set  of  more  detailed  algorithm 
descriptions  in  Section  6.  Six  numerical  examples  are  presented  in  Section  7  to  illustrate  the 
performance  of  the  schemes.  Finally,  Section  8  lists  some  generalizations  and  conclusions. 
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2.2  Relevant  Facts  from  Approximation  Theory 

The  principal  tool  of  this  paper  is  a  somewhat  detailed  analysis  of  Fourier  series  of  functions 
4> :  [-7r,7r]  — >  C  given  by  the  formula 

<t>{x)  =  e-*®'  •  (9) 

where  b>  \  and  c  are  real  numbers.  We  present  this  analysis  in  the  Lemmas  and  Theorems  of 
this  subsection,  numbered  2.5-2.10. 

Lemmas  2.5  and  2.6,  provide  two  inequalities  which  are  used  in  Theorem  2.7.  Theorems  2.7- 
2.9  are  intermediate  results  leading  to  Theorem  2.10,  which  explains  how  to  approximate  func¬ 
tions  of  the  form  using  a  small  number  of  terms,  and  is  the  principal  result  of  this  section. 
For  all  approximations  we  derive  error  bounds  which  allow  us  to  perform  numerical  computa¬ 
tions  to  any  specified  accuracy. 


Lemma  2,5  For  any  real  b  >  ^,c  and  any  integer  k, 

|2  cos((c  -  k)x)dx  +  e-'^  •  <  27re-^  •  ^1  +  .  (10) 

Proof.  Using  the  triangle  inequality  and  Lemma  2.4  we  have 
I2  r  •  cos((c  -  A:)x)dx  +  e*(‘=-^)^dx  < 


< 

too 

2  e- 

Jtt 

< 

2Tre-^^ 

< 

2ire-^ 

Lemma  2.6  For  any  real  b>  \,c  and  any  integer  k, 

^2  J  e~^^^  cos((c  -  A:)x)dx  +  •  J  e'^‘^~^^^dx  <  ■*" 

Proof.  Integrating  by  parts  we  have 

2  /  •  cos((c  -  k)x)dx 

Jtt 

=  — ^  sin((c  -  fc)!)]”  + sin((c  -  A:)i)dx  (13) 

c  —  a:^  c  —  k  Jt 

= - sin((c  —  fc)?r)  H — f  xe~^^^  s\n{{c  —  k)x)dx. 

c  —  k  C  —  k  J-jr 

After  rearranging  the  terms  in  (13)  and  integrating  by  parts  again  we  obtain 


2  [  e  cos{{c  —  k)x)dx  +  — — —  sin((c  — /:)7r) 
Jr  C  —  k 


4b  r°° 

C  —  k  J'jr 

4b 

~{c-ky 
4b  f 
(c  -  k)^  V 
4b  / 
(c  -  ky  I 


xe  sin((c  -  k)x)dx 


^  xe  cos((c  —  fc)x)j  —  J  (1  — 26i^)e  cos{{c  —  k)x)dx'^  (14) 


^:re  +  j  ®  *^^dx  +  J  x-2bxe  ’’^^dx^ 
(ne-’’^"  +  y’°°e-^^*dx+ [-xe-^^']~  + 


e-*’^"dx 
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Finally,  due  to  (14)  and  Lemmas  2.1  and  2.4  we  have 

U  f°°  g-6x2  ^  _  r  ^i(c-k)x^, 

1  A  J— 


46e-^  /  2  ^ 

-  <  (TTW-r  +  ^j 


4be-^ 

(c-fc)2- 

SiTre-**^" 

(c-ky 


(l  +  ^).  (15) 


The  following  theorem  provides  an  explicit  expression  for  the  coefficients  of  a  Fourier  series 
which  approximates  functions  of  the  form  (9). 

Theorem  2.7  Let  4){x)  =  for  any  real  b  >  |,c.  Then,  for  any  x  6  (— tt,  tt), 


4>{x)-  f;  p,e'^^  < 


where 


for  k  =  -00, ...  ,00. 

Proof.  We  denote  by  Ok  the  k-th  Fourier  coefficient  for  (/>,  so  that  for  x  €  (-Jr,  tt), 

4>ix)=  £  (18) 

k=—oo 

and  due  to  Lemma  2.3  and  equation  (17)  we  have 
<^k  =  4>{x)e-'^^dx 

=  pk  -  —  [  cos((c  -  k)x)dx, 

TT  Jv 

Rearranging  equation  (19)  and  applying  Lemmas  2.5  and  2.6  we  obtain  the  inequalities 


riT 

<^k 

-  Pk- 

c 

2-x 

/  €‘^e-’*^dx 

J—ir 

p-bir'^ 

CTk 

-  Pk- 

c 

27r 

/  e‘“e-'*^dx 

/  —  TT 

ic-ky 


4 


and  it  now  follows  from  the  combination  of  (18),  (20)  and  (21)  that,  for  any  x  €  (— tt.tt), 


OO 


ikx  g-fc’T*  .  g»C® 


/:=— OO 
OO 

E  « 

k=—oo 


ikx 


Ok-  Pk~ 


j-6ir* 

2t 


f 


^icxg-ikxdx 


fc,|c-fc|>3 


fc.Jc— A:l<3 


Some  elementary  analysis  yields 

^  ¥^9^  Jz  ■^2-9  +  3“9 


k=3 


and  substituting  (23)  into  (22)  we  have 


k=~oo 


<  e 


(-!)• 


To  complete  the  proof  we  make  use  of  the  triangle  inequality  and  (24)  to  obtain 


OO 

OO 

+  \e-^^  .  e’“^| 

4>(x)-  Y 

< 

•Kx)-  Y 

k^—oo 

fc=— OO 

(22) 


(23) 


(24) 


(25) 

□ 


Remark  2.1  According  to  Theorem  2.7,  functions  of  the  form  e"*’*^*e*“  can  be  approximated 
by  a  Fourier  series  whose  coefficients  are  given  analytically,  and  the  error  of  the  approximation 
decreases  exponentially  as  b  increases. 


Remark  2.2  The  coeflficients  pk  as  defined  by  (17)  have  a  peak  a.t  k  =  [c],  the  nearest  integer 
to  c,  and  decay  exponentially  as  k  ±oo.  We  keep  only  the  g  +  1  largest  coefficients,  where 
the  integer  q  is  chosen  such  that 

q  >  Ab-K,  (26) 


so  as  to  satisfy  the  inequality 

g-(9/2)V46  <  (27) 
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The  following  theorem  estimates  the  truncation  error  under  the  conditions  of  Remark  2.2 
and  thus  provides  a  way  of  approximating  functions  of  the  form  (9)  by  a  g-term  series. 

Theorem  2.8  Let  4>{x)  =  e*®®  for  any  real  b  >  ^,c,  and  let  q  be  an  even  integer  such 

that  q  >  4bx.  Then,  for  any  x  E  (— 7r,7r), 

W+9/2 


*:=W-9/2 

where  {pk}  are  defined  by  (17). 

Proof.  For  any  x  E 


<  .{4b  +  9), 


[c]+9/2 

‘^(®)  ~  E  Pk^''"^ 

< 

+ 

E  pk^''^^ 

+ 

E  pk^'’’^ 

*:=W-9/2 

fc=— oo 

*>W+9/2 

k<lc]-q/2 

Due  to  (17)  and  the  triangle  inequalty  we  have  the  inequalities 

~  ^-(c-kflAb 


fc>lc]+?/2 


<  E 


oo  g-fc^/4i> 


E  Pk^ 

,fc<[cl-q/2 


.ikx 


k=lc]r,/7+i  2v^  "  ^/2^  ’ 

W-9/2-1  g_(c_fc)2/4i  oo  ^-k^/4b 


7«  +  r  .  f  1  +  , 

Jq/2  \  2ql2/ 


fe=-oo  ‘  V 

Some  elementary  analysis  and  an  application  of  Lemma  2.4  yields 

CO 

E»- 

k=q/2 

and  it  follows  from  the  combination  of  (26),  (27)  and  (32)  that 

<.-  .(i  +  i). 

Substituting  (33)  into  (30)  and  (31)  we  have 

2e-^’^* 


v/2?  ■ 


OO 

E' 

k=ql2 


—k^/4b  ^  iff* 


E_  „ikx 

Pke 

k>[c]+q/2 


+ 


E 

k<[c]-q/2 


y/^ 


^  9  ’ 


and  finally,  substituting  (25)  and  (34)  into  (29),  we  obtain 


[c]+?/2 

-  E  Pk^'’'^ 

k=[c]-ql2 


<  e' 


(28) 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 


(35) 

□ 


The  following  corollary  describes  a  scheme  for  approximating  using  a  series  of  q  terms. 


Corollary  2.9  Suppose  that  m  >  2  is  an  integer  and  that  the  conditions  of  Theorem  2.8  are 
satisfied.  Then,  multiplying  both  sides  of  (28)  by  we  obtain 

H+9/2 

g.cr_g6x2.  ^  -(46  +  9) 

fc=[c]-9/2 

<  (36) 

for  anyx€ 

Finally,  Theorem  2.10  makes  use  of  a  simple  linear  scaling  to  generalize  the  inequality  (36) 
from  to  any  interval  [—d,d\. 

Theorem  2.10  Let  b  >  ^,c,d>  0  be  real  numbers,  and  let  m>2,q>  4&7r  be  integers.  Then, 
for  any  x  £  [— <f,  d\, 

[cm(i/ir]+9/2 

gicx  _  g6(ix/m<i)2  .  ^  ^  ^-bx^(l-l/m^)  .  (4J  4.  9)  (37) 

k^[cmd/ir]—q/2 

where  {pk}  are  defined  by  (17). 

Remark  2.3  The  estimated  error  bounds  obtained  in  the  above  theorems  are  rather  pes¬ 
simistic.  Numerical  estimates  for  the  actual  errors  can  be  found  in  Appendix  A  to  this  paper. 

3  Exact  Statement  of  the  Problem 

In  the  remainder  of  this  paper  we  will  operate  under  the  following  eissumptions: 

1 .  a;  =  {wqj  •  •  • »  ^n}  and  x  =  {iq,  ...,xn}  are  finite  sequences  of  real  numbers. 

2.  uk  €  [-iV/2,  Ar/2]  for  k  =  0,...,N. 

3.  Xj  £  [-TT,  7r]  for  j  =  0, . . . ,  N. 

4.  or  =  {ao,.  .  ■,Qn},  f  =  {/-/V/2>--  -  0  =  {0-N/2y  ■  ■  ■>  Pn/2}  1  9  =  {^O,  ■  •  -  .Pn}, 

7  =  {70, . . . ,  7;v}  and  h  =  {/iq,  . . . ,  h/v}  are  finite  sequences  of  complex  numbers. 

We  will  consider  the  problems  of  applying  and  inverting  the  matrix  of  the  Fourier  kernel  and 
its  transpose,  i.e.  we  are  interested  in  the  transformations  F,G  :  and  their 

inverses  defined  by  the  formulae 

N 

/,.  =  F(a),  =  (38) 

fc=o 

N/2 

9:  =  G{0h=  E  (39) 

k=-N/2 
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defined  by  the 
(40) 


•  Problem  1 

Given  a,  find  /  =  F{a). 

•  Problem  2 

Given  /3,  find  g  =  G{P). 

•  Problem  3 
Given  7,  find  h  = 

•  Problem  4 

Given  /,  find  a  = 

•  Problem  5 

Given  g,  find  P  =  G~^{g). 


Remark  3.1  If  Xj  =  — Wj  •  27r/N,  then  G  =  F*. 

We  will  also  consider  the  more  general  transformation  H  : 
formula  ^ 

ik=0 

More  formally,  we  consider  the  following  problems 


Remark  3.2  We  wish  to  perform  aU  calculations  with  a  fixed  relative  accuracy  £  >  0.  In  the 
case  of  Problem  1,  for  instance,  we  are  looking  for  a  vector  /  =  {f-N/2,  •  •  • ,  In/t}  such  that 


II/- /II 

ll/ll 


<  £. 


In  this  sense,  aU  algorithms  described  in  this  paper  are  approximate  ones. 


(41) 


4  Informal  Descriptions  of  the  Algorithms 

4.1  Algorithms  1,  2  and  3  for  Problems  1,  2  and  3 

Observation  4.1  According  to  Theorem  2.10,  any  function  of  the  form  e'“  can  he  accurately 
represented  on  any  finite  interval  on  the  real  line  using  a  small  number  of  terms  oj  the  form 
and  this  number  of  terms,  q,  is  independent  of  the  value  of  c. 

The  FFT  algorithm  applies  the  matrix  of  the  Fourier  kernel  to  arbitrary  complex  vectors 
in  O(fVlogiV)  operations  when  are  integers  and  {ij}  are  equally  spaced  in  [— x,7r].  For 
the  efficient  application  of  the  transformations  described  by  (38),  (39)  and  (40),  we  relate  these 
more  general  cases  to  the  equispaced  case  of  the  FFT.  Observation  4.1  is  used  in  two  ways  to 
achieve  this: 
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•  To  approximate  each  e*"*®  in  terms  of  a  9-term  Fourier  series. 

•  To  approximate  the  value  of  a  Fourier  series  at  each  Xj  in  terms  of  the  values  of  this  series 
at  the  nearest  q  equispaced  nodes. 

This  interpolation  between  equispaced  and  non-equispaced  sets  of  points  can  thus  be  per¬ 
formed  in  0(Nq)  operations.  The  overall  complexity  of  each  of  Algorithms  1,  2  and  3  will 
therefore  be  0(iV  log  iV  Nq)  operations. 

4.2  Algorithms  4  and  5  for  Problems  4  and  5 

Here  we  are  interested  in  applying  the  complex  matrices  A~^  and  to  arbitrary  complex 

vectors  where  the  elements  of  A  are  defined  by 

Ajk  =  (42) 

for  j  =  0,..  .,N  and  k  =  —N/2, . .  .,N/2. 

We  make  use  of  the  following  two  simple  observations. 

Observation  4.2  The  matrix  AA*  is  Toeplitz,  and  furthermore,  its  2N  -b  1  distinct  elements 
can  he  computed  using  Algorithm  1. 

Proof.  It  is  obvious  from  (42)  that 

N  N 

{AA’)jk  =  X!  (43) 

/=o  1=0 

which  is  a  function  only  of  {j  —  k),  and  is  of  the  same  form  as  (38),  the  description  for  Problem 
1.  □ 


Observation  4.3  From  elementary  matrix  identities  we  see  that 

A-^  =  A*(AA*)-^  (44) 

(A*)-^  =  (AA*)-M.  (45) 

The  Toeplitz  matrix  AA*  can  be  applied  to  arbitrary  vectors  in  0{N  log  N)  operations 
using  an  FFT-based  discrete  convolution.  (AA*)“^  can  therefore  be  applied  to  a  vector  in 
0(«:(A)  •  Alog  N)  operations  using  the  conjugate  gradient  method  where  k(A)  is  the  condition 
number  of  A. 

Observation  4.4  As  application  of  A*  and  A  are  0{N  log  N  +  Nq)  procedures  using  Algo¬ 
rithms  1  and  2,  A~^  end  (A*)“^  can  be  applied  in  0{k{A)  •  A  log  A  -f  Nq)  operations  due  to 
Observation  4-3. 

Remark  4.5  It  is  well  known  that  the  condition  number  of  A  is  1  if  the  points  {xj}  are  equally 
spaced.  While  the  condition  number  deteriorates  as  the  distribution  of  points  becomes  more 
non-uniform,  in  many  cases  of  practical  interest  the  points  will  be  fairly  uniformly  spaced,  so 
the  condition  number  will  not  be  very  large. 
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4.3  Algorithm  6  for  a  Variant  of  Problem  5 

The  following  lemma  describes  a  way  of  computing  the  coefficients  of  an  y-term  Fourier  series 
which  is  tabulated  at  N  points. 


Lemma  4.1  Suppose  that  the  N  +  1  function  values  go,---igN  o.re  given  by  the  formula 

N/4 

9i=  E  (46) 

k=-N/4 

and  the  vector^  =  {^o,-  ■  -  Is  the  unique  solution  of  the  linear  system  described  by  the 
equations 


Yi  \  = 

1  0  otherwii 
j=o  1 


otherwise 

fork=  Then,  fork  =  -N  jA,. .  .,N  JA, 

A =i:«i 

1=0 

Proof.  Substituting  for  gj  from  (46)  we  have  for  k  =  —N/A, .  ..,N/A 
N 

3=0 


(47) 


(48) 


N  NI4 

=  E  E 

(49) 

3=0  l=-N/4 

N/4  N 

(50) 

o 

II 

1 

II 

=  Pk. 

(51) 

□ 

Remark  4.6  According  to  (47)  and  Lemma  2.2, 

(62) 

MO  27r  J.^ 

for  k  =  —Nf2,...,N/2.  Thus,  the  set  of  numbers  {^j}  can  be  considered  as  weights  which 
integrate  exactly  all  TV-th  order  trigonometric  polynomials  at  the  nodes  {ij}. 

Observation  4.7  Rewriting  (47)  in  matrix  notation  we  see  that 

=  (0,... ,0,1,0,. ..,0)^  =  A*e  (53) 

where  Ajk  =  so  the  vector^  is  real,  and  can  be  computed  using  the  algorithm  for  Problem 
4  as  introduced  in  Section  4-2. 
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Observation  4.8  Equation  (48)  is  of  the  same  form  as  (38).  Thus,  provided  the  vector  ^  is 
known,  the  vector  /?  can  be  computed  in  0{N  log  N  +  Nq)  operations  using  the  algorithm  for 
Problem  1  as  introduced  in  Section  4.1. 

Remark  4.9  According  to  Observation  4.8,  if  a  function  described  by  an  A/2-term  Fourier 
series  is  tabulated  at  N  arbitrary  nodes,  the  N/2  coefficients  can  be  obtained  in  0{N  log  N+Nq) 
operations. 

Also,  due  to  Observations  4.7  and  4.4,  the  precomputation  of  the  numbers  {^j}  needed  for 
this  algorithm  requires  0(k(A)  •  A  log  A'  +  Nq)  operations. 


5  Notation 


In  this  section  we  introduce  the  notation  to  be  used  in  the  detailed  descriptions  of  the  algorithms 
in  the  next  section. 

For  an  integer  m  >  2  and  a  real  number  6  >  0,  we  will  define  a  real  number  £  >  0  by 

£^g-fc»"(i-i/m").(4j  +  9),  (54) 

and  we  will  denote  by  q  the  smallest  even  natural  number  such  that 

q  >  4bT.  (55) 

For  an  integer  m  and  a  set  of  real  numbers  {w*}  we  will  denote  by  pk  tbe  nearest  integer 
to  muk  for  k  =  0,. .  .,N,  and  by  {Pjk}  a  set  of  real  numbers  defined  by  the  formula 

P-k  =  — - (50) 

2'\/6t 

for  A:  =  0, . . . ,  iV  and  j  =  -q/2, .  ■.,q/2. 


Observation  5.1  Setting  d  =  ir  in  Theorem  2.10  we  see  that 

9/2 

gtu/fcr  _  g6(z/Tri)*  .  ^  p.^  . 


<  £ 


(57) 


1  j=-9/2  1 

for  any  k  =  0, . . .,  N  and  any  x  €  [— tt, t],  where  e  is  defined  by  (54). 

For  a  given  set  of  complex  numbers  we  will  denote  by  {rj}  the  unique  set  of  complex 
coefficients  such  that 


N  9/2 


k=\  3=-qn 


mAf/2-1 


j=—mN/2 


(58) 
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so  that 


(59) 


^  OCk-  Pjk- 

j,k,lik+3=l 

We  will  denote  by  {Tj}  a  set  of  complex  numbers  defined  by  the  formula 

mN/2— 1 

Tj=  ^  .  ^2^ikj/mN  (go) 

k=-mN/2 

for  j  =  — miV/2, . . . ,  mN  12  -  1. 

Further,  we  will  denote  by  {/,}  auother  set  of  complex  numbers  defined  by  the  form\ila 

fj  =  gilS’ri/mA')*  .  JV  (gl) 

for  j  =  -  N/2, . . . ,  N/2. 

Observation  5.2  Combining  (57)  -  (61)  with  the  tr'angle  inequality,  we  see  that 

(62) 

k=0 

for  j  —  —N/2, . .  .,N/2,  where  {fj  =  F(a)j]  are  defined  by  (38). 

For  an  integer  m  and  a  set  of  real  numbers  {a;^}  we  will  denote  by  i/j  the  nearest  integer  to 
XjmN/2ir  for  j  =  0, . . . ,  N,  and  by  {Qjk}  a  set  of  real  numbers  defined  by  the  formula 

Qjk  =  (63) 

2v^ 

for  J  =  0, . . iV  and  k  =  —q/2, . .  .,q/2. 


Observation  5.3  Setting  d  =  iV/2  in  Theoreni  2.10  we  see  that 

7/2 

_  ^h(2-KklmNY  .  ^  Q.^  .  gi[vj+l)2irklmN 

/=-,/2 


<  e 


(64) 


for  any  j  =  0, . .  .,N  and  any  k  €  [— A'^/2, iV/2],  where  e  is  defined  by  (54). 

For  a  given  set  of  complex  numbers  {/3k},  we  will  denote  by  {u^}  a  set  of  complex  numbers 
defined  by  the  formula 

Uk  =  Pk-  (g5) 

for  k  =  —N/2, . . .,  N/2,  and  by  {[/;}  a  set  of  complex  numbers  defined  by  the  formula 

N/2 

U{=  (66) 

k=-N/2 
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for  I  =  -miV/2, . . . ,  mN  12  -  1. 

Further,  we  will  denote  by  {gj}  another  set  of  complex  numbers  defined  by  the  formula 


l=-q/2 


for  j  =  0,...,N. 


Observation  5.4  Combining  (64)  -  (67)  with  the  triangle  inequality,  we  see  that 

N 

l5j-5il  (68) 

k=0 

for  j  =  Q,..  .,N ,  where  {gj  =  G{0)j)  are  defined  by  (39). 

For  a  set  of  real  numbers  {xj}  we  will  denote  by  rjj  the  nearest  integer  to  XjN /2t:  for 
j  =  0, . . IV,  and  by  a  set  of  real  numbers  defined  by  the  formula 

R  _  .  g-(i,iV/2^-(r,j+*))V46  (09) 

2v6x 

for  j  =  0, . . . ,  fV  and  k  =  -q/2, . . . ,  q/2. 

Observation  5.5  Setting  d  =  N/2  in  Theorem  2.10  we  see  that 

9/2 

gikxjim  _  ^b{2irk/mN)^  .  ^  r.^  .  gi(rij+l)2irk/mN  ^  ^  (JQ) 

J=-,/2 

for  any  j  =  0,. .  .,N  and  any  k  6  [— iV/2, iV/2]  where  e  is  defined  by  (54). 

For  a  given  set  of  complex  numbers  {7*},  we  will  denote  by  {uj}  the  unique  set  of  complex 
coefficients  such  that 

N  9/2  mN/2 

E7fc-  E  E  Vj-e'^-^^,  (71) 

k=0  j=-g/2  j=—mN/2 

SO  that 

VI  =  E 

},k,Vk+i=l 

We  denote  by  {F;}  a  set  of  complex  numbers  defined  by  the  formula 


k=-mN/2 
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for  /  =  —m^N/2, . . . ,  w?N/2  —  1. 

Further,  we  will  denote  by  {hj}  another  set  of  complex  numbers  defined  by  the  formula 

9/2 

E  RirVr,,^i  (74) 

'=-9/2 

for  j  =  0, . .  .,iV. 

Observation  5.6  Combining  (57)  and  (70)  -  (74)  with  the  triangle  inequality,  we  see  that 

*  N 

|/i,-A,l</).El7Jtl  (75) 

fc=o 

for  j  =  0, . .  .,N ,  where  {hj  =  -ff (7);}  are  defined  by  (40),  and 

6  =  26-*”^' .  (46  +  9).  (76) 

For  a  set  of  real  numbers  {xj},  A  will  denote  a  complex  matrix  whose  elements  are  given 

by 

Ajk  =  (77) 

for  k  =  -N/2, .  ..,N/2  and  j  =  0,...,N,  and  a  =  {o_;y,. .  .,aiv}  will  denote  a  set  of  complex 
numbers  defined  by  the  formula 

=  (78) 

j=0 

Finally,  ^  =  {^oj  •  •  • ,  Civ}  will  denote  a  real  vector  defined  by 

(yl*)-i(0,.. .,0,1,0,. ..,0)^.  (79) 

Remark  5.7  It  is  clear  from  Observation  4.2  that 

{AA^)jk  =  fli-fc.  (80) 
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6  Detailed  Descriptions  of  the  Algorithms 

This  section  contains  step  by  step  descriptions  and  operation  counts  for  the  six  algorithms  of 
this  paper.  In  the  tables  below  we  will  make  use  of  the  facts  that  q  ~  log(^)  and  -C  N. 

Algorithm  1. 

Step  Complexity  Description 

Init  0{Nq)  Comment  [Input  parameters  are  the  vector  {wo» •  • aJid  a 

real  number  e  >  0.] 

Choose  b  and  q  in  terms  of  £ 
do  k  =  Q,N 

Determine  /ijt,  the  nearest  integer  to  muk 
do  j  =  -?/2,?/2 

Compute  Pjk  according  to  (56) 
end  do 
end  do 

do  j  =  -NJ2,N/2 
Compute 
end  do 

1  0(Nq)  Comment  [Input  parameter  is  the  vector  {oq,  . . . ,  ojv}-] 

Comment  [Compute  Fourier  coefficients  rj.] 

do  /:  =  0,  JV 

do  j  =  -q/2,qf2 

^  '^Uk+J  "b 

end  do 
end  do 

2  0{rnN\o^N)  Comment  [Evaluate  this  Fourier  Series  at  equispaced  points  in 

[— m7r,jn7r]  using  inverse  FFT.] 

Compute  Tj  =  Zk=IlN/2  for  j  =  -mN/2, . .  .,mN/2  - 

3  Comment  [Scale  the  values  at  those  points  which  lie  in  [-T,  tt].] 

do  j  =  ~NI2,NJ2 
fj  =  Tj  ■ 
end  do 

Total  0(7V’ ■  log(i)  +  miV  •  log  iV) 
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Algorithm  2. 

Step  Complexity 


Description 


Init  0(Nq)  Comment  [Input  parameters  are  the  vector  {xo5-*-i*a?}  and  a 

real  number  e  >  0.] 

Choose  6  and  q  in  terms  of  e 
do  j  =  0,N 

Determine  Uj,  the  nearest  integer  to  XjmNf2'K 
do  A:  =  -9/2, 9/2 

Compute  Qjk  according  to  (63) 
end  do 
end  do 

do  fc  =  -NI2,NI2 
Compute 
end  do 

1  0{N)  Comment  [Input  parameter  is  the  vector  {j3_a?/2»-  • 

Comment  [Compute  new,  scaled  Fourier  coefficients.) 

do  k  =  -N/2,N/2 
Uk  =  Pk’ 
end  do 

2  0(mN\ogN)  Comment  [Evaluate  this  Fourier  Series  at  equispaced  points  in 

[-7r,7r]  using  inverse  FFT.j 

Compute  Uj  =  for  j  =  —mN  12, . .  .,mlV/2  —  1. 

3  0{Nq)  Comment  [Compute  approximate  values  at  desired  points  in 

terms  of  the  values  at  equispaced  points.) 

do  j  =  Q,N 

do  k  =  -9/2,972 

^  Pj  +  Qjk  ■ 

end  do 
end  do 

Total  OimN  ■]ogN  +  N  ■  log(i)) 
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Algorithm  3. 

Step  Complexity  Description 

Init  0{Nq)  Comment  (Input  parameters  are  the  vectors  {uq,  . .  .,u}n}  and 

{iQ, . .  and  a  real  number  £  >  0.] 

Choose  b  and  q  in  terms  of  £ 
do  k  =  0,N 

Determine  //jt,  the  nearest  integer  to  muk 
do  3  =  -9/2, 9/2 

Compute  Pjk  according  to  (56) 
end  do 
end  do 

do  A:  =  —mN/2,mNl2 
Compute 
end  do 
do  j  =  OjN 

Determine  rjj,  the  nearest  integer  to  XjN/2n 
do  A:  =  —972,9/2 

Compute  Rjk  according  to  (69) 
end  do 
Compute 
end  do 

1  0{Nq)  Comment  [Input  parameter  is  the  vector  {70, . . . ,  7n}-] 

Comment  (Compute  Fourier  coefficients  v,  .] 
do  A:  =  0,A 

do  j  =  -9/2, 9/2 

^Vifc+i  ^  'Ik 

end  do 

end  do 

2  O(mN)  Comment  [Scale  the  coefficients.] 

do  A:  =  —mN/2,mN/2 
Vki-Vk- 
end  do 

3  0{m‘^N  log  N)  Comment  (Evaluate  this  Fourier  Series  at  equispaced  points  in 

(— m7r,m;r]  using  inverse  FFT.j 

Compute  Vj  =  Y^=imNi2  for  j  —  -m'^N/2, . . . ,  Tn?N/2  - 
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4  0(Nq)  Comment  [Compute  approximate  values  at  desired  points  in 

terms  of  the  values  at  equispaced  points.] 
do  j  =  0,N 

do  A:  =  —q/2,ql2 

hj  ■*—  hj  +  Rjk  ■  Vrij+k 

end  do 
end  do 

5  0(N)  Comment  [Scale  the  values.) 

do  j  =  0,N 

hj  ^  hj  ■ 
end  do 

Total  0{m?N  -logJV  +  iV  •log(i)) 

Algorithm  4. 

Step  Complexity  Description 

Init  0(NlogN  +  Nq)  Initialization  for  Algorithm  1. 

Initialization  for  Algorithm  2. 

Compute  elements  {ajt}  of  Toeplitz  matrix  as  defined  by 

(78)  using  Algorithm  1. 

1  0(N  log  N  +  Nq)  Compute  A/  using  Algorithm  2. 

2  0(k(A)  •  N  log  A)  Compute  d  =  (AA*)”^(A /)  using  Conjugate  Gradient  algorithm. 

Total  0(k{A)  ■  N -logN  +  N -logil)) 

Algorithm  5. 

Step  Complexity  Description 

Init  0{N\ogN  +  Nq)  Initialization  for  Algorithm  1. 

Compute  elements  {a*}  of  Toeplitz  matrix  (AA*)”^  as  defined  by 
(78)  using  Algorithm  1. 

1  0(/c(A)  •  N  log  A)  Compute  {AA‘)~^g  using  Conjugate  Gradient  algorithm. 

2  0{N\ogN  +  Nq)  Compute  =  A*((AA*)“^5)  using  Algorithm  1. 

Total  0(k(A) -A -logA  +  iV  •log(l)) 
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Algorithm  6. 

Step  Complexity 

Description 

Init 

0{k{A)-NIosN  +  Nq) 

Initialization  for  Algorithm  5. 

Compute  ^  as  defined  by  (79)  using  Algorithm  5. 

1 

0{N) 

Compute  gj  =  ^jgj  for  j  =  0,. . N . 

2 

OiN  log  N  +  Nq) 

Compute  P  =  A'g  using  Algorithm  1. 

Total 

OiN  -logN  +  N  -logil)) 

The  storage  requirements  of  an  algorithm  are  also  an  important  characteristic.  From  the 
above  descriptions  for  the  initialization  steps  the  asymptotic  storage  requirements  for  each 
algorithm  are  of  the  form 

X-N-q  (81) 

where  the  coefficient  A  is  softwaxe-  and  hardware-dependent. 

7  Numerical  Results 

VVe  have  written  FORTRAN  implementations  of  the  six  algorithms  of  this  paper  and  have  tested 
these  programs  on  the  Sun  SPARCstation  1  for  a  variety  of  input  data.  Six  such  examples  are 
presented  in  this  section,  one  for  each  algorithm. 

Tables  1-6  contain  accuracies  and  CPU  timings  for  the  algorithms  with  computations  per¬ 
formed  in  both  single  and  double  precision  arithmetic,  and  the  input  size  N  varying  between 
64  and  4096.  In  addition,  each  table  contains  the  CPU  times  required  to  solve  the  same  set  of 
problems  via  a  direct  calculation,  and  Tables  1-3  include  timings  for  an  FFT  of  the  same  size. 
Tables  1-3  also  contain  the  accuracies  of  the  direct  single  precision  calculations. 

Two  measures  of  accuracy  were  chosen  for  each  example.  In  Examples  1,  2  and  3,  these  are 
defined  by  the  formulae 


and 


max 

0<j<N 


\h  -  f:\ 


E2 


N 


Ni=o 


(82) 


(83) 


where  a  is  the  input  vector,  /  is  the  result  of  a  direct  computation  in  double  precision  arithmetic 
and  /  is  the  result  of  the  computation  being  considered. 
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In  Examples  4,  5  and  6,  they  are  given  by 


and 


where  a  is  the  input  for  a  direct  double  precision  computation,  and  a  is  the  result  of  applying 
the  algorithm  to  the  result  of  this  computation. 


Remark  7.1  The  formulae  (82)  -  (85)  measure  fairly  accurately  the  errors  of  all  single  preci¬ 
sion  computations.  However,  they  can  only  provide  rough  estimates  of  the  errors  produced  by 
the  double  precision  versions  of  the  algorithms. 

Several  technical  details  of  our  implementations  appear  to  be  worth  mentioning  here: 

1.  Each  implementation  consists  of  two  main  subroutines:  the  first  is  an  initialization  stage 
in  which  the  matrix  operators  of  the  sJgorithm  are  precomputed  and  stored,  and  the 
second  is  an  evaluation  stage  in  which  these  operators  are  applied.  Successive  application 
of  the  linear  transformations  to  multiple  vectors  requires  the  initialization  to  be  performed 
only  once. 

2.  For  the  single  precision  versions  of  each  algorithm  our  choice  of  parameters  was  m  =  2, 
h  —  0.5993  and  q  =  10.  For  double  precision  we  chose  m  =  2,b  =  1.5629  and  q  =  28. 

3.  The  algorithms  as  described  in  this  paper  all  require  an  FFT  of  size  proportional  to  N, 
and  thus  will  perform  efficiently  whenever  the  FFT  does.  This  restriction  on  the  FFT  size 
can  be  removed  by  means  of  a  simple  modification,  thereby  ensuring  that  the  algorithms 
perform  efficiently  for  any  choice  of  N.  In  our  implementations,  these  changes  were  made. 

4.  In  the  direct  methods  for  Problems  1  and  2,  we  used  the  fact  that  =  (e'^-’)*'  to 
reduce  the  number  of  exponential  computations  from  iV^  to  N. 

5.  iV^  exponentials  are  required  in  the  direct  method  for  Problem  3,  and  for  larger  N, 
the  available  memory  on  the  machine  is  insufficient  for  the  precomputation  and  storage 
of  these  numbers.  The  direct  implementation  we  used  for  this  problem  computes  each 
exponential  when  it  is  needed. 

6.  Standard  LINPACK  Gaussian  Elimination  subroutines  were  used  as  the  direct  methods 
for  comparing  timings  in  Examples  4,  5  and  6.  Estimated  timings  are  presented  for  larger 
N,  where  this  computation  became  impractical. 
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Example  1, 

Here  we  consider  the  transformation  F  :  of  Problem  1  as  defined  by  the  formula 

=  (86) 

fc=0 

for  j  =  -N/2, . .  .,N/2.  In  this  example,  {wq.  •  • -.wjv}  were  randomly  distributed  on  the 
interval  [—N/2,  N/2]  and  {qq,  .  •  .jOtv}  were  generated  randomly  on  the  unit  square  in  the 
complex  plane  defined  by  the  formulae 

0  <  Re(z)  <  1,  0  <  Im(2)  <  1.  (87) 

The  results  of  applying  Algorithm  1  to  this  problem  are  presented  in  Tables  1(a)  and  1(b). 


Table  1(a) 

Example  1,  Single  Precision  Computations 


N 

Errors 

Timings  (sec.) 

Algorithm 

Direct 

Algorithm 

Direct 

EFT 

P 

E2 

E2 

Init. 

Eval. 

64 

0.175  E-05 

0.190  E-05 

0.337  E-06 

0.811  E-06 

0.011 

0.003 

KEM 

128 

0.973  E-06 

0.203  E-05 

0.449  E-06 

0.162  E-05 

0.022 

0.006 

mm 

^  1 

256 

0.113  E-05 

0.348  E-05 

0.100  E-05 

0.374  E-05 

0.044 

0.014 

m 

512 

0.138  E-05 

0.602  E-05 

0.168  E-05 

0.769  E-05 

0.088 

0.030 

0.49 

1  £ 

1024 

0.280  E-05 

0.128  E-04 

0.193  E-05 

0.133  E-04 

0.174 

0.067 

1.90 

m 1 

2048 

0.267  E-05 

0.243  E-04 

0.424  E-05 

0.268  E-04 

0.348 

0.141 

7.47 

m  1 1 

4096 

0.551  E-05 

0.453  E-04 

0.532  E-05 

0.565  E-04 

0.701 

0.323 

30.24 

Table  1(b) 

Example  1,  Double  Precision  Computations 


■1 

Errors 

Timings  (sec.) 

■i 

Eca 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

FFT 

64 

0.495  E-14 

0.634  E-14 

0.036 

0.02 

0.001 

128 

0.689  E-14 

0.104  E-13 

0.075 

WBm 

0.06 

0.002 

256 

0.717  E-14 

0.119  E-13 

0.148 

■EH 

0.22 

0.005 

512 

0.306  E-14 

0.164  E-13 

0.297 

0.84 

0.012 

1024 

0.460  E-14 

0.310  E-13 

0.600 

3.23 

0.026 

2048 

0.694  E-14 

0.625  E-13 

1.204 

■«■ 

12.55 

0.059 

4096 

0.129  E-13 

0.126  E-12 

2.418 

49.69 

0.132 
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Example  2. 

Here  we  consider  the  transformation  G  :  ^  of  Problem  2  eis  defined  by  the  formula 

N/2 

GW)j=  £  (88) 

k=-N/2 

for  j  =  0, N .  In  this  example,  were  randomly  distributed  on  the  inter’val 

[— 7r,7r]  and  {P-n/2t  •  •  Pn/t}  were  generated  randomly  on  the  unit  square  in  the  complex 
plane  defined  by  the  formulae 


0  <  Re(2)  <  1,  0  <  lm(0)  <  1. 


(89) 


The  results  of  applying  Algorithm  2  to  this  problem  are  presented  in  Tables  2(a)  and  2(b). 


Table  2(a) 

Example  2,  Single  Precision  Computations 


N 

Errors 

Timings  (sec.) 

Algorithm 

Direct 

Algorithm 

Direct 

EFT 

Eoo 

E2 

Eoo 

E2 

Init. 

Eval. 

64 

0.121  E-05 

0.173  E-05 

0.237  E-06 

0.506  E-06 

0.011 

0.003 

1S9 

0.0008 

128 

0.595  E-06 

0.337  E-05 

0.362  E-06 

0.985  E-06 

0.022 

0.007 

lEl 

0.0015 

256 

0.133  E-05 

0.620  E-05 

0.233  E-06 

0.158  E-05 

0.043 

0.015 

WSi 

0.0034 

512 

0.953  E-06 

0.890  E-05 

0.828  E-06 

0.334  E-05 

0.086 

0.033 

0.0078 

1024 

0.164  E-05 

0.996  E-05 

0.431  E-05 

0.523  E-05 

0.174 

0.067 

0.0174 

2048 

0.223  E-05 

0.156  E-04 

0.687  E-05 

0.824  E-05 

0.345 

0.149 

7.11 

0.0352 

4096 

0.470  E-05 

0.338  E-04 

0.759  E-05 

0.134  E-04 

0.687 

0.332 

29.03 

0.0846 

Table  2(b) 

Example  2,  Double  Precision  Compuations 


N 

Errors 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

EFT 

64 

0.249  E-14 

0.814  E-14 

0.038 

■BB 

0.01 

128 

0.501  E-14 

0.746  E-14 

0.075 

Wmm 

0.05 

256 

0.418  E-14 

0.623  E-14 

0.150 

0.18 

512 

0.356  E-14 

0.831  E-14 

0.297 

0.69 

1024 

0.793  E-14 

0.192  E-13 

C.598 

2.72 

2048 

0.138  E-13 

0.405  E-13 

1.192 

BB 

10.86 

ifiiaSI 

4096 

0.278  E-13 

0.904  E-13 

2.417 

43.36 

0.132 
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Example  3. 

Here  we  consider  the  transformation  H  :  of  Problem  3  as  defined  by  the  formula 

N 

=  (90) 

k=Q 

foT  j  —  0, N .  In  this  example,  {wq,  . .  were  randomly  distributed  on  the  interval 

[—N/2,  N/2],  {xq,  . . . ,  lyv}  were  randomly  distributed  on  the  interval  [— tt,  tt]  and  {70, . . . ,  7^^} 
were  generated  randomly  on  the  unit  square  in  the  complex  plane  defined  by  the  formulae 

0  <  Re(z)  <  1,  0  <  Im(2)  <  1.  (91) 

The  results  of  applying  Algorithm  3  to  this  problem  are  presented  in  Tables  3(a)  and  3(b). 


Table  3(a) 

Example  3,  Single  Precision  Computations 


N 

Errors 

Timings  (sec.) 

Algorithm 

Direct 

Algorithm 

Direct 

EFT 

Eqo 

E2 

Eca 

E2 

Init. 

Eval. 

64 

0.111  E-05 

0.291  E-05 

0.570  E-06 

0.707  £-06 

0.022 

0.007 

0.24 

|IKIIIIj» 

128 

0.204  E-05 

0.405  E-05 

0.593  E-06 

0.145  E-05 

0.044 

0.015 

0.63 

256 

0.127  E-05 

0.450  E-05 

0.900  E-06 

0.296  E-05 

0.087 

0.032 

2.60 

0.0034 

512 

0.181  E-05 

0.683  E-05 

0.111  E-05 

0.490  E-05 

0.177 

0.067 

10.63 

0.0078 

1024 

0.304  E-05 

0.176  E-04 

0.221  E-05 

0.126  E-04 

0.352 

0.145 

43.60 

0.0174 

2048 

0.482  E-05 

0.356  E-04 

0.260  E-05 

0.264  E-04 

0.708 

0.327 

176.75 

0.0352 

4096 

0.734  E-05 

0.751  E-04 

0.389  E-05 

0.530  E-04 

1.404 

0.702 

714.40 

0.0846 

Table  3(b) 

Example  3,  Double  Precision  Computations 


N 

Errors 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

FFT 

64 

0.222  E-13 

0.320  E-13 

0.078 

0.21 

1 

128 

0.247  E-13 

0.370  E-13 

0.153 

0.85 

256 

0.249  E-13 

0.334  E-13 

0.305 

KB 

3.43 

1 

512 

0.145  E-13 

0.232  E-13 

0.605 

13.62 

HR  1 

1024 

0.237  E-13 

0.416  E-13 

1.210 

54.71 

HR  1 

2048 

0.194  E-13 

0.795  E-13 

2.403 

WSBM 

219.38 

HR' 

4096 

0.411  E-13 

0.120  E-12 

4.827 

1.588 

889.04 

0.132 
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Example  4. 

Here  we  consider  Problem  4  of  Section  3.  In  this  example,  the  numbers  were  defined  by 
the  formula 

Wfc  =  +  (fc  +  0.5  +  6k)  ■  ^  ^  (92) 

for  /:  =  0, . .  .,iV,  where  Sk  were  randomly  distributed  on  the  interval  [—0.1,0.!].  In  addition, 
the  numbers  {qq,  . .  .,o:Ar}  were  generated  randomly  on  the  unit  square  in  the  complex  plane 
defined  by  the  formulae 

0  <  Re(2)  <  1,  0  <  Im(^)  <  1,  (93) 

and  the  numbers  {/_Ar/2)  •  •  ■ifN/2}  were  computed  directly  in  double  precision  arithmetic  ac¬ 
cording  to  the  formula 

=  (94) 

k=o 

The  vector  /  wa.s  then  used  as  input  for  Algorithm  4.  Results  of  this  experiment  are  presented 
in  Tables  4(a)  and  4(b). 


Table  4(a) 

Example  4,  Single  Precision  Computations 


N 

Errors 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

64 

0.492  E-05 

0.339  E-05 

0.02 

0.05 

0.36 

128 

0.154  E-04 

0.568  E-05 

0.04 

0.11 

2.78 

256 

0.424  E-04 

0.128  E-04 

0.08 

0.20 

23.0 

512 

0.809  E-04 

0.252  E-04 

0.18 

0.45 

184 

1024 

0,203  E-03 

0.474  E-04 

0.36 

0.82 

1472  (est.) 

2048 

0.428  E-03 

0.979  E-04 

0.78 

1.91 

11776  (est.) 

4096 

0.106  E-02 

0.195  E-03 

1.66 

4.64 

94208  (est.) 

Table  4(b) 

Example  4,  Double  Precision  Computations 


N 

Errors 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

64 

0.143  E-13 

0.109  E-13 

0.07 

0.17 

0.37 

128 

0.208  E-13 

0.149  E-13 

0.11 

0.34 

2.96 

256 

0.493  E-13 

0.256  E-13 

0.20 

0.75 

23.6 

512 

0.121  E-12 

0.500  E-13 

0.48 

1.67 

189 

1024 

0.279  E-12 

0.926  E-13 

0.94 

3.55 

1512  (est.) 

2048 

0.593  E-12 

0.192  E-12 

1.95 

7.75 

12096  (est.) 

4096 

0.138  E-11 

0.375  E-12 

4.02 

18.28 

96768  (est.) 
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Example  5. 

Here  we  consider  Problem  5  of  Section  3.  In  this  example,  the  numbers  {xj}  were  defined  by 
the  formula 


— X  +  2r  • 


j  +  0.5  +  Sj 
N  +  1 


(95) 


for  j  =  0, . .  .,iV,  where  Sj  were  randomly  distributed  on  the  interval  [— 0.1,0.!].  In  addition, 
the  numbers  {P-n/2i-  ■  -  iPn/t}  ''^®re  generated  randomly  on  the  unit  square  in  the  complex 
plane  defined  by  the  formulae 


0  <  Re(2)  <  1,  0  <  Im(2)  <  1,  (96) 

and  the  numbers  {go, . . .  ,g;v}  were  computed  directly  in  double  precision  arithmetic  according 
to  the  formula 

N/2 

9j=  (97) 

h=-N/2 

The  vector  g  was  then  used  as  input  for  Algorithm  5.  Results  of  this  experiment  are  presented 
in  Tables  5(a)  and  5(b). 


Table  5(a) 

Example  5,  Single  P.cCioion  Computations 


N 

Errc.'s 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

64 

0.212  E  05 

0183  E-05 

0.02 

0.06 

128 

0.117  E-04 

0.523  E-05 

0.04 

0.10 

2.78 

256 

0.190  E-04 

0.955  E-05 

0.09 

0.19 

23.0 

512 

0.287  E-04 

0.202  E-04 

0.19 

0.41 

184 

1024 

0.560  E-04 

0.376  E-04 

0.37 

0.83 

1472  (est.) 

2048 

0.106  E-03 

0.752  E-04 

0.78 

1.90 

11776  (est.) 

4096 

0.225  E-03 

0.148  E-03 

1.66 

4.67 

94208  (est.) 

Table  5(b) 

Example  5,  Double  Precision  Computations 


N 

Errors 

Timings  (sec.) 

E<x> 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

64 

0.310  E-13 

0.120  E-13 

0.06 

0.18 

0.37 

128 

0.389  E-13 

0.146  E-14 

0.10 

0.36 

2.96 

256 

0.577  E-13 

0.204  E-13 

0.24 

0.76 

23.6 

512 

0.673  E-13 

0.325  E-13 

0.47 

1.61 

189 

1024 

0.118  E-12 

0.817  E-13 

0.97 

3.54 

1512  (est.) 

2048 

0.19C  E-12 

0.134  E-12 

1.86 

7.73 

12096  (est.) 

4096 

0.429  E-12 

0.288  E-12 

3.93 

18.21 

96768  (est.) 
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Example  6. 

Here  we  consider  the  variant  of  Problem  5  which  was  described  in  Section  4.3.  In  this  example, 
the  numbers  {ij}  were  defined  by  the  formula 


Xj  =  —TT  +  27r 


j  +  0.5  +  6j 
N  +  1 


(98) 


for  j  =  0,...,N,  where  6;  were  randomly  distributed  on  the  interval  [—0.1,0.!].  In  addition, 
the  numbers  {i9_A74,  •  • . ,  A7V/4}  were  generated  randomly  on  the  unit  square  in  the  complex 
plane  defined  by  the  formulae 


0  <  Re(2)  <  1,  0  <  Im(2)  <  1,  (99) 

and  the  numbers  {go?  •  •  •  »yiv}  were  computed  directly  in  double  precision  arithmetic  according 
to  the  formula 

N/4 

9:=  (100) 

k=-N/4 

The  vector  g  was  then  used  as  input  for  Algorithm  6.  Results  of  this  experiment  are  presented 
in  Tables  6(a)  and  6(b). 


Table  6(a) 

Example  6,  Single  Precision  Computations 


N 

Errors 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

64 

0.195  E-05 

0.174  E-05 

0.09 

0.004 

0.36 

128 

0.412  E-05 

0.286  E-05 

0.17 

0.007 

2.78 

256 

0.105  E-04 

0.913  E-05 

0.34 

0.014 

23.0 

512 

0.195  E-04 

0.147  E-04 

0.70 

0.031 

184 

1024 

0.474  E-04 

0.320  E-04 

1.41 

0.066 

1472  (est.) 

2048 

0.836  E-04 

0.631  E-04 

3.05 

0.147 

11776  (est.) 

4096 

0.173  E-03 

0.135  E-03 

7.22 

0.330 

94208  (est.) 

Table  6(b) 

Example  6,  Double  Precision  Computations 


N 

Errors 

Timings  (sec.) 

-  -o 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

64 

0.878  E-14 

0.762  E-14 

0.32 

0.008 

0.37 

128 

0.102  E-13 

0.981  E-14 

0.62 

0.018 

2.96 

256 

0.213  E-13 

0.180  E-13 

1.25 

0.039 

23.6 

512 

0.444  E-13 

0.376  E-13 

2.64 

0.083 

189 

1024 

0.679  E-13 

0.520  E-13 

5.82 

0.170 

1512  (est.) 

2048 

0.151  E-12 

0.115  E-12 

12.37 

0.359 

12096  (est.) 

4096 

0.308  E-12 

0.232  E-12 

27.61 

0.766 

96768  (est.) 
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The  following  observations  can  be  made  from  Tables  1-6  above,  and  are  in  agreement  with 
results  of  our  more  extensive  experiments. 

1.  The  errors  produced  by  Algorithms  1,  2  and  3  are  comparable  with  those  produced  by 
the  corresponding  direct  methods. 

2.  The  timings  for  Algorithms  1  and  2  are  similar,  which  is  to  be  expected  since  Problem  2 
is  the  adjoint  of  Problem  1.  Algorithm  3  is  about  twice  as  costly,  which  is  in  agreement 
with  the  fact  that  it  is  a  synthesis  of  Algorithms  1  and  2. 

3.  In  single  precision  computations  for  this  particulcir  architecture,  implementation  and 
range  of  A,  Algorithms  1  and  2  are  roughly  4  times  as  costly  as  an  FFT  of  the  same  size. 
For  double  precision  the  ratio  is  roughly  6.  These  ratios  decrease  as  N  increases. 

4.  The  extrapolated  break-even  point  of  Algorithms  1  and  2  is  at  roughly  A  =  32  if  the 
initialization  time  is  ignored.  If  the  initialization  time  is  included,  the  break-even  point  is 
at  A  =  256.  For  Algorithm  3,  the  break-even  points  are  at  A  =  8  without  initialization 
and  at  A  =  32  with  initialization. 

5.  The  timings  for  Algorithms  4  and  5  are  similar  as  expected,  since  Problem  5  is  the  adjoint 
of  Problem  4. 

6.  The  break-even  points  of  Algorithms  4  and  5  are  at  roughly  A  =  32.  For  Algorithm  6, 
the  break-even  points  are  at  A  =  32  if  the  initialization  time  is  taken  into  account,  and 
at  A  =  16  if  it  is  ignored. 

7.  Algorithms  1  and  2  tend  to  be  slightly  more  accurate  than  their  inverses,  Algorithms  4 
and  5. 

8.  The  initialization  for  Algorithm  6  is  computationally  costly,  but  subsequent  evaluations 
require  much  less  CPU  time  than  evaluations  for  Algorithm  5. 

Remark  7.2  The  CPU  timings  for  Algorithms  1,  2  and  3  are  independent  of  the  particular 
distributions  of  u?  and  x,  whereas  the  timings  for  Algorithms  4  and  5  are  sensitive  to  the 
distributions  of  these  vectors. 
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8  Generalizations  and  Conclusions 

The  results  of  this  paper  can  be  generalized  in  the  following  ways: 

1.  Simple  modifications  to  Algorithms  1,  2  and  3  will  allow  the  efficient  application  of  the 
linear  transformations  -♦  defined  by 


N 


fc=0 

N/2 

for  j  =  -M/2,...,Af/2, 

(101) 

Gmi 

=  iZ  Pk- 

k=-N/2 

N 

for  j  =  0, . . .,  M, 

(102) 

k=0 

for  j  =  0, . . . ,  M. 

(103) 

These  changes  have  been  implemented. 

2.  The  algorithms  of  this  paper  also  assume  that  Uk  £  [— iV/2,  N/2]  and  Xj  G  [-^r,  tt].  Other 
distributions  can  be  handled  by  appropriately  partitioning  the  vectors  cj  and  x,  treating 
each  partition  separately  and  finally  combining  the  results.  The  following  observation 
describes  translation  operators  which  can  be  used  for  each  partition,  in  combination  with 
one  of  Algorithms  1,  2  or  3. 

Observation  8.1  Leta,b>  Q,c,d>  0  be  real  numbers  and  suppose  that  Uk  €  [a— 6,o  +  6] 
for  k  —  0,. .  .,N  and  Xj  £  [c  —  d,c-\-  d\  for  j  =  0, . . . ,  Af .  Then  we  can  write 


N 

N 

(104) 

jt=0 

fc=0 

fc=0 

(105) 

where 

a'k  = 

^'k  =  (wfc-a)d/7r,  (106) 

x'j  =  (xj  _c)7r/d€  [-7r,7r]. 

Remark  8.2  Such  an  algorithm  will  perform  efficiently  when  the  points  within  a  parti¬ 
tion  are  close  together  and  there  are  very  few  partitions,  and  not  so  efficiently  if  the  points 
are  widely  separated  and  there  are  many  partitions.  Most  cases  likely  to  be  encountered 
in  practice  fall  in  the  former  category. 
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3.  A  paper  describing  a  set  of  algorithms  based  on  a  different  interpolation  technique  is 
currently  in  preparation. 

4.  One  of  the  more  far-reaching  extensions  of  the  results  of  this  paper  is  a  set  of  algorithms 
for  problems  in  higher  dimensions.  Investigations  into  this  are  currently  in  progress. 

5.  The  Helmholtz  equation  in  2  dimensions  is  given  by 

(107) 

and  hzis  particular  solutions  of  the  form 

=  (108) 

where  fP  +  =  k^.  Solutions  of  this  equation  consist  of  linear  combinations  of  such 

functions  subject  to  some  boundary  conditions,  and  the  results  of  this  paper  admit  a 
generalization  which  constitutes  a  fast  Helmholtz  solver. 

In  conclusion,  a  group  of  algorithms  has  been  presented  for  the  rapid  application  and 
inversion  of  matrices  of  the  Fourier  kcinel.  These  problems  can  be  viewed  as  generalizations 
of  the  Discrete  Fourier  Transform,  and  the  algorithms,  while  making  use  of  certain  simple 
results  from  analysis,  are  very  versatile,  and  have  wide-ranging  potential  applications  in  many 
branches  of  mathematics,  science  and  engineering. 


Appendix  A 


In  this  section  we  present  numerical  estimates  of  the  errors  in  approximating  using  q  terms 
of  the  form  for  the  values  of  m,  b  and  q  used  by  the  algorithms.  As  q  is  independent 

of  the  choice  of  c,  we  only  considered  the  case  c  =  0.  The  results  are  presented  in  Table  7. 
The  approximation 

?/2 

k=-q/2 

to  the  constant  function  4i(x)  =  1  with  pk  defined  by  equation  (17)  was  computed  at  n  =  1000 
equally  spaced  nodes  Xk  in  and  the  entries  in  the  table  are  defined  as  follows: 


•  Eoo  is  the  maximum  absolute  error  defined  by  the  formula 


=  max  l<^(xfc)  -  (i>ixk)\. 

l<k<n 


•  E2  is  the  relative  L2  error  defined  by  the  formula 


E2  = 


(110) 


(111) 


•  jEj3  is  the  error  bound  of  Theorem  2.10  defined  by  the  formula 

=  +  (112) 


Table  7 


m 

6 

9 

E<x> 

E2 

Eb 

2 

0.5993 

10 

0.825  E-05 

0.176  E-05 

0.135  E-00 

2 

1.5629 

28 

0.400  E-13 

0.580  E-14 

0.163  E-03 

30 
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